Embedding on to a one-dimensional crystal 

J.E. Inglesfield 
School of Physics and Astronomy, Cardiff University, 
The Parade, Cardiff, CF24 3AA, United Kingdom 
(Dated: March 12, 2010) 

Abstract 

A simple expression is derived for the band structure of a one-dimensional periodic potential in 
terms of two solutions of the Schrodinger equation within the unit cell, one with a zero-derivative 
boundary condition on the left-hand end of the cell and the other with zero derivative on the right- 
hand end. Prom this starting point, a new expression is derived for the embedding potential - this 
can be added to the Hamiltonian for the surface region of a crystal to replace the semi-infinite 
substrate, in a one-dimensional approximation. The results are demonstrated in calculations of the 
band structure and embedding potential for Al in the [001] direction, and the surface electronic 
structure of the Al(OOl) surface. 

PACS numbers: 71.15.-m, 71.20.-b, 73.20.-r 



I. INTRODUCTION 



Embedding provides a way of including the effects of the substrate in a calculation of 
electronic structure over a restricted region of space. ^ For example, embedding allows us 
to find the electronic structure at a surface by solving the Schrodinger equation in the 
surface region, adding embedding potentials to the Hamiltonian to include the effects of 
the semi-infinite substrate and vacuum regions.'^ The embedding potentials ensure that 
the wave-functions (or the Green function) in the surface region have the correct boundary 
corrections on the boundaries of the region, without matching wave-functions explicitly. The 
alternative to embedding is to solve the Schrodinger equation for a slab of material, or a 
periodic array of slabs. However, slab calculations do not give the energy continuum of 
bulk states, and unless the slab is very thick, the localised surface states interact across the 
slab. Embedding calculations have neither of these drawbacks. Several methods have been 
developed over the years for finding the embedding potential for a semi-infinite substrate, 
with the full crystal potential, for solving the energy-dependent Schrodinger equation in full- 
potential calculations of surface electronic structure.'^ In this paper we develop a very fast 
method for calculating the embedding potential to replace a semi-infinite one-dimensional 
periodic potential. 

Going from the full three-dimensional crystal potential to a one-dimensional potential may 
seem a retrograde step, but Chulkov et a/.'^have recently developed one-dimensional models 
of the bulk and surface potential (figure 1), which give an accurate description of electronic 
states at the surface. This has proved particularly useful in many-body and lifetime studies 
of Shockley and image potential-induced surface states.l^ This one- dimensional potential is 
of course felt by electrons moving in three dimensions, with a wave-function of the form 

^(r) = exp(iK.R)V(^), (1) 

where K is the free-electron wave-vector parallel to the surface, which lies in the R 
plane. Our own interest lies in using these model surface potentials in the time- dependent 
Schrodinger equation, which we solve by a new embedding technique.'^ The crystal substrate 
is replaced by a time-dependent embedding potential, essentially the Fourier transform of the 
energy-dependent embedding potential. This must be evaluated at a very fine energy grid 
over a very wide energy range - hence we need a method for finding the energy- dependent 
embedding potential as efficiently as possible. 
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FIG. 1. One-dimensional potential for modelling the Al(OOl) surface.lSl At Zy the surface region is 
embedded on to the vacuum Coulomb tail, and at Zc it is embedded on to the semi-infinite crystal 
potential. 

The general idea of embedding is that we partition space into two (or more) regions, 
solving the Schrodinger equation explicitly in what we call region I, with the rest of space, 
region II, replaced by an embedding potential added to the Hamiltonian of region I. The 
embedding potential is defined over the interface S between regions I and 11^^ Using a 
variational method, it can be shown that the wave-function for the system, ip{r), satisfies 
the following Schrodinger equation in region I, 



__V2 + ^^(r) ) ^(r) + 5{r - 
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where the integral is over S!^ The embedding potential S is a function of two interface 
coordinates and is evaluated at trial energy e; the energy derivative term gives E at the 
required energy E, to first order in [E — e). The embedding potential gives the generalised 
logarithmic derivative over S of the wave-function in region 11,1^ 
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and this ensures that ijj in region I matches correctly in amphtude and derivative on to 
the solution in II. In three-dimensional applications, we find S from the Green function 
for region II satisfying the zero-derivative boundary condition on 5'.'^ However, in the one- 
dimensional application described in this paper, we use Eq. |3]to find the embedding potential 
- in one-dimension, this becomes a relation between amplitude and derivative at the point 
between regions I and II. 

As the use of S suggests, the embedding potential is a form of self-energy, an energy- 
dependent, possibly complex potential added on to the Hamiltonian to replace a region of 
phase space or Hilbert space. In fact the embedding potential in the linear combination 
of atomic orbitals formalismP' is usually called the self-energy, often in connection with 
conduction through linear molecules attached at each end to metallic contacts.'^^IEl xhe con- 
tacts and associated electron reservoirs are replaced in the tight-binding Hamiltonian of the 
system by self-energies (embedding potentials in tight-binding guise In our embedding 
method, in which space is partitioned, we can use any convenient basis set for expanding 
the wave-function or Green function in region I. 

We now briefiy describe the structure of this paper. In section II we shall use solutions 
of the Schrodinger equation in one unit cell to calculate the complex band structure, which 
describes the allowed solutions for the semi-infinite crystal at energy E, and the logarithmic 
derivative of these solutions, hence the embedding potential S(-E). Related results have been 
derived by Butti,'''- using a different method, and his expression for the one-dimensional em- 
bedding potential was used in a recent analysis of electron spectroscopies from adsorbates.'^ 
The result for the band structure has previously been derived by Kohn,^ - and with particular 
reference to the sinusoidal potential (the Mathieu problem), by McLachlan.^ The expres- 
sions are remarkably simple, and the band structure formula in particular may be useful 
in teaching, where one-dimensional potentials such as the Kronig-Penney model frequently 
serve as an introduction to band theory.l^^l In section HI we shall illustrate the use of the 
embedding potential in a calculation of the density of states at the Al(OOl) surface using 
the one-dimensional model potential of Chulkov et al!^ The expression for the embedding 
potential of a one-dimensional crystal is useful not only in the energy domain, but we are 
also using it in Fourier transform (with an extra factor of 1/ to study time- dependent 
processes at surfaces, such as electron emission. 

Atomic units are used in this paper, with e"^ = h = rUe = I- 
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II. BAND STRUCTURE AND EMBEDDING POTENTIAL 



Our starting point for finding tlie band structure and embedding potential is the Green 
function formula for the wave-function ip{z) in some interval in terms of the derivatives at 
the ends of the range, 

i;{z) = ^-[G{z,a)^'{a)-G{z,0)^'m, < z < a, (4) 

where we take the range to be the unit cell between z = and z = a of the infinite one- 
dimensional crystal; G is the Green function with zero-derivative boundary conditions at the 
ends of the unit cell. This formula is analogous to the equation in electrostatics giving the 
potential inside some region of space in terms of the boundary values of the electric field,'^!] 
and can be derived in exactly the same way. 

This equation satisfied by the wave-function within the unit cell does not depend on tp{z) 
in the rest of the crystal, but we now impose the Bloch form of wave-function. Evaluating 
ip a.t z = 0, Eq. |4] becomes 

iP{0) = ^tp\0)[G{0,a)exp{ika) - G{0,0)], (5) 



where k is the Bloch wave- vector. But we can also evaluate Eq. |4] at 



z = a, ffivmE 



ip{a) = -ilj'{a)[G{a,a) — G{a,0) exp{—ika)]. (6) 

Now the logarithmic derivative ip'/ip is invariant to a lattice displacement, so comparing 
these two equations we obtain 

^(0, a) exp{ika) - G(0, 0) = G{a, a) - G{a, 0) exp{-ika). (7) 

The Green function is symmetric in its spatial variables, so G{0,a) = G{a,0), and Eq. [T] 
simplifies to 

G{0,0) + G{a,a) 
= 2G(0,a) • 

This can be simplified further by using the expression for the Green function in terms of 
wave-functions (f)i{z) and 02 (-z), which satisfy the boundary conditions of the Green function 
at each end of the unit cell: (pi satisfies the Schrodinger equation with (f)'i{0) = 0, and at the 
other end of the unit cell 02(a) = 0, 

py(0i,02) 



where W is the Wronskian.'^ Furthermore, we take 0i(O) = 1, 02 (o) = 1. Then substituting 
Eq. |9] into Eq. [8] gives the remarkably simple expression for the wave- vector, 

, . 0i(a) + 02(O) . . 

cos{ka) = . (10) 

So we can find the Bloch wave-vector corresponding to a particular energy by two integra- 
tions through the unit cell, one in each direction, to find (pi and 02- If the potential is 
symmetric with respect to the origin, 0i(a) = 02(0), and the result simplifies to 

cos(A;a) = 0i(a), (11) 

Using wave-function matching, Kohri^S! gives a result related to Eq. 10, though instead 



of using 02; his second wave-function has zero amplitude at z = and unit derivative. 
Abramowitz and StegurP^ quote Eq. 11 in section 20.3.10 in the chapter on Mathieu func- 



tions, and this is proved by McLachlan!!^ using a Fourier expansion of the wave- function. An 
alternative approach is given by ButtiP he uses the same starting wave-functions as Kohn 
to find the transfer matrix; this relates solutions of the Schrodinger equation at each end of 
the unit cell, and its eigenvalues give the Bloch phase factors. 

If the solutions 0i and 02 at real energy E satisfy |0i(a) + 02(O)|/2 > 1, then E lies in a 



bulk band-gap, and the wave-vector k satisfying Eq. 10 is complex (either pure imaginary, 
or complex with a real part at the Brillouin zone boundary). This corresponds to a Bloch 
solution forbidden in the infinite bulk crystal, but allowed in the semi-infinite crystal, for 
example in the case of a crystal with a surface.'^ More generally we can find the band- 



structure at complex E - we shall need this in section III. In this case, k satisfying Eq. 10 
is always complex, and the two solutions correspond to two different physical cases. One 
solution corresponds to waves travelling to the right and decaying in this direction, with 
Imfc > 0; the other solution, with ImA; < is travelling and decaying to the left. Which 
solution we take depends on the particular problem - in our case, with the semi-infinite 
bulk crystal lying in the positive z-direction (figure [T]) the physical solution corresponds to 
Im > 0. 

To illustrate the method, we evaluate the complex band structure of a one-dimensional 
pseudopotential corresponding to Al in the [001] direction, 

V{z) = Acos{27rz/a), (12) 
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FIG. 2. Complex band structure for Al in the [001] direction with the one-dimensional potential. 
ReE is plotted along the y-axis, and ImE is kept fixed at 0.002 a.u. Solid line, Rek; dashed line, 
Im k. 

where the amphtude is given by A = 0.0618 a.u. and the lattice constant a is 3.8 a.u. - the 
bulk potential of figure [l|^ This is the Mathieu problem, but our method holds for arbitrary 
potential and arbitrary origin. We calculate the band structure at complex energy, keeping 
the imaginary part of the energy fixed at 0.002 a.u. At each energy we find (pi and 02 by 
integrating the Schrodinger equation through the unit cell using Numerov's method,!^ with 
a spatial interval of 0.002 a.u., starting off the integrations using the method described by 
Quiroz Gonzales and Thompson, The resulting band structure, smoothed by working with 
the small imaginary energy, is shown in figure [2j Re /c is plotted in the extended zone scheme 
to give a continuous curve, weaving its way from band to band around the bulk band gaps; 
as we expect, Imk is non-zero below the bottom of the band, and in the band gaps. Our 
analysis of the complex band structure is important now that we turn to the embedding 
potential. 

We need the embedding potential Ec to replace the semi-infinite bulk crystal, which we 
assume lies to the right of Zc (figure [T]) , where we now put the origin of the unit cell. In the 
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one-dimensional case, Eq. [3] simplifies to the logarithmic derivative, 



(13) 



2 ^(0) ' 

where ip is the solution of the Schrodinger equation for the semi-infinite crystal travelling 
or decaying as z ^ oo^ We can find this logarithmic derivative directly from Eq. [s} and 
substituting the Wronskian expression for the Green function, Eq. |9j we obtain the following 
result for the embedding potential, 

2[exp(zfca) -02(0)]' ^ ^ 

The wave-vector k in this expression corresponds to the wave travelling or decaying to the 
right - this is precisely the value of k which was discussed above in the band structure 
calculation (figure [2|. This result is different in form from the expression derived by Butti,'^ 
but the two must be identical. 
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FIG. 3. Embedding potential for a bulk Al substrate with the one-dimensional potential, evaluated 
at Zc = 10 a.u., as a function of Re-E. ImE is kept fixed at 0.0002 a.u. Solid line, ReSc; dashed 
line, ImSc. 

We now calculate the embedding potential for the Al substrate at Zc = 10 a.u., at a 
complex energy with an imaginary part equal to 0.0002 a.u. The only additional feature 
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compared with the band structure calculation is that the expression for the embedding 
potential involves the Wronskian, in which we evaluate the derivatives by finite differences. 
The results are shown in figure |3| and we note that Im E is negative, as we require from 
causality. 



III. SURFACE DENSITY OF STATES 



The Al(OOl) surface is suitable for demonstrating this embedding potential - the one- 
dimensional crystal potential of Chulkov et al.^ should work well in this s-p bonded metal, 
and the (001) surface shows a Shockley surface state,'^ as well as structure induced by the 
image potential.'^Sl 

We shall calculate the local density of states e) in the surface region, the charge 
density of states with energy e. This is given by the sum over states 

a{z,e) = Y,\Uz)?^{^-^i). (15) 

i 

where il)i{z) is a wave-function of the system with energy e,, and it can be written in terms 
of the Green function evaluated at an energy with a small imaginary part, 

a{z,e) = —\vsiG{z,z]e + irf). (16) 

TT 

(Note that we use e to denote a real energy, and E to denote an energy which may be 
complex - here, E = e + if].) 

Let us take the surface region, where we calculate G, between z^ on the vacuum side and 
Zc on the crystal side (figure [T]), and then in this region G satisfies the Schrodinger equation 
embedded on both sides, 



Id^G 



+ {V-E)G{z,z';E) + 6{z-Zc) 



2 dz^ 

IdG 



1 3G 

-— + S,(E)G(z„z';E) 



+5{z - Zy 



+E,{E)G{z.,,z';E) 



S{z-z'). (17) 



The surface region is embedded at z^ on to the crystal, with the crystal embedding potential 
Sc, and at z^ on to the vacuum region, which is replaced by S^. Unlike Eq. |2} the embedded 
Schrodinger equation for the Green function does not contain the energy derivative of the 
embedding potentials; this is because we know the energy at which they should be evaluated 
- the energy E of the Green function. 
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FIG. 4. Embedding potential for the vacuum region, evaluated at Zy = —10 a.u. with the image 
plane at zq = —3.44 a.u., as a function of Rei?. \m.E is kept fixed at 0.0002 a.u. Solid line, ReS„; 
dashed line, ImS^,. 

Although it is not the main topic of this paper, we shall say a few words about the 
embedding potential S^,, which replaces the vacuum region outside the crystal. Here the 
electron feels the Coulomb tail of the image potential, 

where Vq is the vacuum zero (we take the average potential in the bulk crystal as the zero of 
energy), and is the position of the image plane. Again we calculate the embedding poten- 
tial from the logarithmic derivative of the outgoing or decaying solution of the Schrodinger 
equation - this is a combination of the regular and irregular Coulomb function d^" * ^^ ' Fq and 
Go (in the notation of Abramowitz and StegurP^)^ with angular momentum L = 



'^{^ = Hq {t], p) = Go{t], p) - iMv, P) 
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(19) 



with arguments given by 

p = ^2{E-V,){z - z,) 
^ A^2{E-Vo)' 

Thompson and BarnetlP^ give a rapidly converging continued fraction expression for 
Hq' /Hq, from which we can immediately find the vacuum region embedding potential 
T,^. Figure |4] shows S^,, at Zy = —10 a.u. with the image plane at zq = —3.44 a.u. as a 
function of energy, taking ImE = 0.0002 a.u. The structure just below the vacuum zero at 
E = 0.577 a.u. comes from bound states of the Coulomb potential. 

We are now in a position to solve the embedded Schrodinger equation in region I, the 
near-surface region, using a basis set expansion for G, 

G{z,z';E) = Y,G,,{E)x^{z)x,{z'). (21) 

The Schrodinger equation Eq. 17 then reduces to a matrix equation, 

{Hij + T,ij{E) — ESij) Gjk = Sik, (22) 

j 

where the Hamiltonian matrix is given by 

H^J = I r dz^^ + r dzx.{z)V{z)xUz), (23) 
the embedding matrix by 

S,,(E) = S,(E)x.(^e)x,(^c) + Y.,{E)x^{z,)xA^.), (24) 
and S is the overlap matrix. 



Sij = / dzxi{z)xj{z). (25) 

J Zv 

We usually use trigonometric basis functions, 

.cos m even ^ ^ 

Xm{^)={ '''' , (26) 

sin m odd 

where ( is measured from the mid-point of region I, 

C = z-^ (27) 
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FIG. 5. Surface density of states of Al(OOl), the local density of states integrated through region 
I, taken between Zy = —10 a.u. and = +10 a.u. 



and a value of D > [z^ — z^) /2 gives a range of logarithmic derivatives at z^ and Zc for 



matching on to the embedding potentials. Solving Eq. [22] then gives us G, and via Eq. [15 
the local density of states. 

Integrating the local density of states over the surface region, we obtain the surface 
density of states for Al(OOl), shown in figure [s] Here we take region I between z^ = —10 a.u. 
and Zc = +10 a.u., and 20 basis functions are used, defined with D = 13 a.u. These results, 
which correspond to K = in three dimensions, show the Shockley surface state near the 
bottom of the band gap, broadened by the imaginary part of the energy at which the the 
Green function and embedding potentials are evaluated, here taken to be 0.0002 a.u. Less 
familiar is the structure just below the vacuum edge, at E = 0.577 a.u., which comes from 
surface resonances induced by the image potential.'^ These are the image potential surface 
states, broadened into resonances by interacting with the continuum of bulk states. All the 
features in the surface electronic structure - the continuum of bulk states, the Shockley and 
image surface states - depend on a correct treatment of the substrate and vacuum regions, 
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via the embedding potentials. 



IV. CONCLUDING REMARKS 

The motivation for this work was to derive the embedding potential for a one-dimensional 
model of a crystal as efficiently as possible, for subsequent use in energy-dependent or time- 
dependent surface calculations. We have shown how the band structure and embedding 
potential can be found by integrating through the unit cell at a chosen energy, starting at 
each end of the cell with a zero-derivative boundary condition. Using these results, we have 
calculated the band structure and embedding potential for Al in the [001] direction as a 
function of complex energy, and applied this to calculating the surface density of states at 
the Al(OOl) surface. 

We are now using these results in studies of time-dependent excitation processes at sur- 
faces, taking the Fourier transform of the embedding potential to find the time-dependent 
embedding potential. But one-dimensional periodic potentials occur in a variety of contexts 
- not least in teaching - and we hope that the results and their derivations will be of wider 
interest. 
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